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I.  INTRODUCTION 


The  high  power  requirements  of  future  spacecraft  missions  present  many 
technological  challenges  to  the  designers  of  spacecraft  thermal  management  systems 
(SCTMS).  The  need  to  buffer  thermal  loads,  such  as  during  burst  power  conditions, 
is  presently  a  major  problem.  One  possible  solution  is  the  application  of  latent  heat 
storage  with  its  use  of  phase  change  materials  (PCM)  to  buffer  the  large  heat  input 
conditions.  The  present  investigation  was  motivated  by  the  need  for  predictive  tools 
necessary  for  the  design  of  SCTMS  utilizing  latent  thermal  energy  storage.  In  this 
research,  a  latent  heat  storage  unit  is  configured  as  shown  in  Fig.  1.  The  PCM  is 
contained  in  a  long  horizontal  tube  assembled  within  a  cylindrical  shell.  The 
temperature  controlled  heat  transfer  fluid  (HTF)  is  pumped  through  the  shell-and-tube 
heat  exchanger  to  release  or  absorb  thermal  energy  to  or  from  the  PCM. 


Figure  1.  Latent  Heat  Storage  Unit  and  Coordinate  System 


Basically,  the  analysis  of  heat  transfer  in  this  latent  heat  storage  unit  is  divided 
into  two  problems:  the  melting/ freezing  problem  of  the  PCM  in  a  horizontal  tube  and 
the  convective  heat  transfer  problem  of  the  HTF  in  an  annular  passage. 

For  the  present  melting/freezing  problem,  the  unknown  tube  wall  temperature  is 
assumed  to  be  time  and  location  dependent.  This  problem  is  rather  difficult  to  solve 
due  to  the  moving  boundary  and  the  unknown  wall  temperature.  By  introducing  the 
enthalpy  method,  this  problem  becomes  much  simpler.  The  governing  equation  of  the 
enthalpy  method  is  similar  to  the  single  phase  conduction  equation.  Using  the 
enthalpy  method  to  approach  a  phase  change  problem  has  the  following  advantages: 
(1)  there  are  no  conditions  to  be  satisfied  at  the  phase  change  boundary;  (2)  there  is 
no  need  to  accurately  track  the  phase  change  boundary;  (3)  there  is  no  need  to  consider 
the  solid  and  the  liquid  sides  separately;  and  (4)  it  allows  a  mushy  zone  between  solid 
and  liquid  phases.  The  unknown  tube  wall  temperature  between  the  heat  transfer  fluid 
and  the  PCM  is  obtained  from  an  energy  balance.  However,  an  iterative  procedure  is 
required.  The  Newton- Raphson  method  is  used  for  this  "root-finding"  problem. 

One  of  the  objectives  of  this  research  is  to  develop  a  numerical  method  to  find 
the  heat  transfer  between  the  HTF  and  the  PCM.  The  comprehensive  understanding 
of  the  heat  transfer  in  the  latent  heat  storage  system  can  provide  the  information  for 
the  design  of  spacecraft  thermal  management  systems.  However,  there  are  many  other 
important  applications,  such  as  metal  casting,  crystal  growth,  and  thermal  energy 
storage  for  building  space  conditioning,  where  the  fundamental  understanding  of  the 
phase  change  phenomena  is  important. 
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II.  BACKGROUND 


A.  MELTING/EREEZING  HEAT  TRANSFER  IN  HORIZONTAL  CYLINDERS 


The  heat  transfer  processes  which  occur  during  melting/ freezing  in  horizontal 
cylinders  have  recently  received  considerable  attention  due  to  a  wide  range  of 
engineering  applications.  Mathematically,  the  problem  belongs  to  the  "moving 
boundary"  class  of  problems,  which  are  difficult  to  find  analytical  solutions.  Webb  and 
co-w'orkers  [I]  studied  the  melting  phenomenology  of  water  in  a  horizontal  cylindrical 
capsule  by  experiment.  The  solid-liquid  interface  and  the  interaction  of  fluid  flow,  due 
to  natural  convection  and  the  motion  of  the  solid  induced  by  the  density  inversion  of 
the  water-ice  system,  were  determined  visually.  Heat  transfer  at  the  solid-liquid 
interface  during  melting  from  a  horizontal  cylindrical  heat  source  embedded  in  a  PCM 
was  studied  by  Bathelt  and  Viskanta  [2].  Two  paraffins,  n-heptadecane,  Tf~  22.2°C 
and  n-octadecane,  Tf  =  28.2 °C  were  used  as  the  PCMs.  The  shape  of  the  solid-liquid 
interface  was  determined  photographically,  and  the  local  heat  transfer  coefficients  were 
measured  using  a  shadowgraph  technique.  Ho  [3]  and  Pannu  [4]  used  the 
vorticity-strcam  function  method  to  deal  with  the  problem  of  phase  change  in 
horizontal  cylinders.  This  method  is  based  on  the  consideration  that  the  heat  transfer 
in  the  liquid  PCM  is  by  natural  convection.  A  perturbation  method  was  used  by  Yao 
[5]  to  solve  the  melting  problem  in  a  horizontal  cylinder.  This  method  is  complex, 
tedious  and  limited  to  small  Stefan  numbers.  Rieger  et  al.  [6,7]  used  a  numerically 
generated  coordinate  system,  which  conformed  to  the  shape  of  the  boundaries,  to  study 
the  melting  process  around  a  horizontal  cylinder.  Due  to  the  changing  sfipc  of  the 
solid-liquid  interface,  this  method  appears  to  be  laborious. 

All  the  aforementioned  studies  did  not  consider  the  temperature  variation  of  the 
PCM  in  the  axial  direction.  Also,  most  are  limited  to  uniform  tube  wall  temperature. 
Experimental  studies  have  also  focussed  on  the  control  of  the  tube  wall  temperature 
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to  be  uniform.  McCabe  [8]  obtained  numerical  solutions  for  time  dependent 
temperature  distributions.  Me  used  the  enthalpy  method  with  an  implicit  scheme  to 
solve  the  melting/freezing  problem  in  cylinders.  However,  only  radial  solutions  were 
obtained.  Like  the  other  authors,  he  also  did  not  take  the  axial  dependence  into 
account.  In  terms  of  the  conjugate  heat  transfer  problems  problems  inmolving  the 
HTF  and  the  PCM  with  "natural"  boundary  conditions  at  the  interface,  several  papers 
[9—11]  have  addressed  the  problem  of  freezing  on  the  outside  of  horizontal  cylinder. 

B.  HEAT  TRANSFER  IN  ANNULAR  PaSSAGLS 

The  annulus  represents  a  common  geometry  employed  in  a  variety  of  heat 
transfer  systems  ranging  from  simple  heat  exchangers  to  the  most  complicated  nuclear 
reactors.  Theoretical  and  numerical  studies  of  heat  transfer  in  an  annulus  can  be 
roughly  divided  into  two  types:  heat  transfer  in  a  hydrodynamically  fully  developed 
region  and  heat  transfer  in  an  entry  region. 

1.  Heat  Transfer  in  Hydrodynamically  Fully  Developed  Region.  For  the  case  of 
an  axially  symmetric  hydrodynamically  fully  developed  fluid  flow  with  constant 
properties  in  a  concentric  annulus,  Reynolds  et  al.  [  12  ]  concluded  that  there  are  four 
types  of  fundamental  boundary  conditions  of  interest:  (1)  a  step  temperature  on  one 
wall  and  the  other  wall  maintained  at  the  inlet  temperature,  (2)  uniform  heat  flux  at 
one  wall  and  the  other  wall  being  kept  adiabatic,  (3)  uniform  wall  temperature  at  one 
wall  and  the  other  wall  remaining  insulated,  and  (4)  uniform  heat  flux  from  one  wall 
with  the  opposite  wall  maintained  at  the  inlet  temperature.  A  very  thorough  treatment 
of  the  annulus  problem  has  been  published  by  Lundberg  et  al.  [  13]  based  on  these 
four  types  of  fundamental  solutions.  They  used  an  eigenvalue  method  to  yield  an  exact 
solution  for  the  entire  thermal  entrance  region,  but  near  the  point  of  a  step-change  in 
cither  the  wall  temperature  or  the  heat  flux  at  the  wall,  a  large  number  of  terms  in  the 


eigenfunction  expansion  are  required.  Worsoe-Schmidt  [  14  ]  obtained  these  four  types 
of  fundamental  solutions  by  using  a  similarity  transformation  method.  The  similarity 
variable  used  was  obtained  from  Leveque  solution.  Hatton  and  Quarmby  [15]  also 
used  the  eigenvalue  method  to  find  solutions  of  heat  transfer  in  annuli.  However,  only 
the  cases  with  cither  uniform  wall  temperature  or  uniform  heat  flux  on  the  inside  wall 
and  with  the  outside  wall  insulated  were  analyzed.  The  case  of  an  arbitrarily  prescribed 
heat  flux  around  the  periphery  of  either  wall,  or  both  walls,  was  studied  by  Sutherland 
et  al.  [  16].  The  technique  used  was  an  expansion  of  the  known  peripheral  heat  flux 
distribution  in  a  Fourier  series.  Because  the  resulting  temperature  distribution  can  also 
be  expressed  as  a  Fourier  series,  the  energy  equation  becomes  a  set  of  ordinary 
differential  equations  in  terms  of  the  eigenfunction.  These  eigenfunctions  were  used 
with  the  Fourier  coefficients  of  any  peripheral  heat  flux  distribution  to  obtain  the 
resulting  wall  temperature  distribution. 

2.  Heat  Transfer  in  Entry  Region.  The  solution  of  the  heat  transfer  problem  with 
simultaneously  development  of  velocity  and  temperature  fields  has  been  the  subject  of 
a  great  deal  of  attention  during  the  past  three  decades.  Due  to  the  fact  that  the 
temperature  profile  cannot  be  obtained  without  knowing  the  velocity  profile,  the 
solution  is  usually  obtained  by  first  solving  the  hydrodynamic  entry  length  problem, 
and  then  solving  the  thermal  entry  length  problem. 

Murakawa  [17]  has  obtained  an  approximate  entrance  length  solution  for  the 
annulus  via  a  series  solution.  The  final  result  only  partially  satisfies  the  boundary 
condition-.;  Lircn  ct  al.  [  18  ]  also  obtained  solutions  for  flow  in  an  annular  region  by 
representing  thf  momentum  equation  by  the  stream  function  and  then  by  using  an 
eigenvalue  n.;'.'  i.  However,  the  solution  is  restricted  to  vanishingly  small  Reynolds 
numbers.  Heaton  al.  [  19  ]  obtained  a  velocity  profile  defined  over  the  entire  region 
by  linearizing  the  momentum  equation  and  solving  the  resulting  equations  analytically. 
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For  the  thermal  entry  length  problem,  Murakawa  [17]  has  obtained  a  solution 
for  an  arbitrary  wall  temperature  at  the  inner  wall  of  the  annulus  by  an  integral 
method.  The  developing  velocity  and  temperature  profiles  used  were  simplified  and 
was  only  applicable  near  the  inner  wall.  The  resulting  solution  did  not  show  all  the 
geometric  effects  nor  did  it  approach  the  fully  developed  solution.  Heaton  et  al. 
[19]  linearized  the  energy  equation  to  find  the  temperature  profile.  Their  solution 
was  obtained  by  an  integral  method.  I  lowcver,  the  only  solution  obtained  was  the  case 
of  uniform  heat  flux  at  one  wall  and  the  other  wall  insulated. 
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III.  ANALYSIS 


A.  ENTHALPY  MODEL  FOR  PHASE  CHANGE  PROBLEM 

For  an  arbitrary  control  volume  V  which  is  fixed  in  space  and  if  conduction  is  the 
only  heat  transfer  mechanism  to  be  considered,  the  time  rate  of  change  of  internal 
energy  has  to  be  equal  to  the  net  rate  at  which  heat  is  conducted  through  V  if  there  is 
no  energy  source  inside  V  and  no  external  work  toward  V.  Mathematically,  this  can 
be  expressed  by  the  following: 


_d_ 

dt 


■f  pUdV  =  f 
Jv  JA 


kVT  .  ndA 


(1) 


Note  that  in  the  absence  of  motion,  the  pressure,  P,  is  independent  of  time.  Thus 


d_ 

dt 


i 


P  dV  =  0 


Since 


pU  =  pi  —  P 

Equation  (1)  can  be  written  in  the  form 


(2) 


(3) 


k  VT 


ndA 


(4) 


Equation  (4)  is  the  general  form  of  enthalpy  equation.  The  i  -  versus  -  T  relationship 
for  the  phase  change  material  is  used  in  conjuction  with  Eq.  (4).  Frequently,  the 
dimensionless  enthalpy,  0,  is  defined  as 
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rv 


6  = 


pV 


p{>  -  ‘t ) 


dV 


(5) 


For  the  case  of  mclting/frcczing  process  in  a  horizontal  tube,  9,  for  the  control 
volume  as  shown  in  Fig.  2  ,  is 


0  = 


P  (i'Jm 


-r'f)Ax 


-P 

'  Jn 


(/-//) 


dV 


(6) 


if  the  buoyancy  effect  in  the  liquid  PCM  is  neglected  and  the  density  is  assumed  to  be 
constant. 

Introducing  9,  Fo,  and  4>,  where 


and 


(7) 


4>  = 


Cp  {T  —  T]) 
k 


(8) 


into  Eq.  (4),  one  has 
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Forward  difference  in  time  and  backward  difference  in  space,  the  explicit 
finite-difference  representation  of  Eq.  (9)  is 


flm+l  —  0<n 
AFo 


Wj+x+Wj-i -*<!>?} 

Ax1 


+  ■ 


'i+ 1 


+  'i 


H4>Zxj-+u)rM-U>Z 


A  r 


(10) 


The  relationships  between  0  and  0  are 


0  =  0 

when 

0  >0 

(11a) 

o 

II 

when 

-  1 <0<O 

(110) 

0  =  0+1 

when 

0  <  —  1 

(11c) 

Equation  (10)  should  be  used  in  conjuctvon  with  the  energy  equation  and  the 
momentum  equation  for  the  HTF,  since  the  conditions  between  the  HTF  and  the  PCM 
are  functions  of  time  and  location.  Due  to  the  unknown  temperature  boundary 
condition  between  the  PCM  and  the  HTF,  the  implicit  scheme  of  the  finite-difference 
equation  for  the  enthalpy  model  appears  to  be  inappropriate.  The  explicit  scheme  is 
therefore  used  since  it  is  the  more  straightforward  method. 

B.  HEAT  TRANSFER  IN  AN  ANNULAR  PASSAGE 

The  current  problem  of  interest  is  heat  transfer  in  an  annular  passage  with  a 
transient  temperature  distribution.  The  dimensionless  differential  equations  associated 
with  this  problem  are  the  continuity  equation,  the  momentum  equation,  and  the  energy 
equation  as  given  below: 
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CONTROL  VOLUME 


Figure  2.  Control  Volume  in  Cylindrical  coordinates 


du  1  d(rv) 

dx  r  dr 


(12) 


du 


1  dp 

2  dx 


+  2(  + 

dr2 


J_  du  . 

r  dr  } 


(13) 


dI.4.XL(ulL+  dT  \  -  1  dT 

dr  2  dx  dr  gr 2  r  dr 


(14) 


These  equations  assume  steady,  laminar,  incompressible  flow  with  constant  fluid 


properties  and  negligible  axial  conduction,  buoyancy  effect,  internal  energy  generation, 


viscous  energy  dissipation,  and  axial  rate  of  change  of  the  radial  shear  stress.  However, 


It' 


for  the  current  problem,  the  flowing  fluid  temperature  is  time  dependent.  Therefore, 
the  the  energy  equation  is  transient.  The  boundary  conditions  for  these  equations  arc 
given  by  the  following: 


u{x,  r2)  =  v(x,  r2)  =  0 

(15) 

u(x,  1)  -  v{x,  1)  =  0 

(16) 

u(0,  r)  =  ue(r) 

(17) 

T(r,  x,  r2)  =  T2( t,  x) 

(  18) 

dT(r,x,  1) 
dr 


(19) 


T(x,0,r)=Te(t) 


(20) 


Equation  (19)  indicates  that  the  shell  is  insulated  at  r—  1.  In  addition  to  these 
geometric  boundary  conditions,  an  initial  boundary  condition  is  required  for  Eq.  (14). 


1.  Hydrodynamically  Fully  Developed  Flow.  If  the  flow  is  hydrodynamically 
fully  developed,  then  the  velocity  is  independent  of  x  and  Eq.  (13)  becomes 


d2 u  1  du 
dr2  r  dr 


1  dp 
4  dx 


constant 


(21) 


The  integration  of  Eq.  (21)  with  the  boundary  conditions  defined  by  Eqs.  (15)  and  (16) 
leads  to  Lamb's  [  20  ]  fully  developed  velocity  profile  for  the  annulus  as  given  below: 
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(22) 


—  r2  +  B  In 


r) 


where 


('2  ~  1) 
In  r2 


and 


M  =  1  +  r\  -  B. 

2.  Flow  in  The  Entrance  Region.  The  problem  of  laminar  flow  heat  transfer  in 
the  entrance  region  was  extensively  studied  during  the  past  three  decades.  However, 
these  studies  were  limited  to  some  simple  cases  such  as  uniform  heat  flux,  uniform  wall 
temperature,  and  steady-state  temperature  distribution.  For  the  current  problem  of  a 
transient  temperature  distribution,  none  of  the  previous  studies  can  be  applied  directly 
and  therefore  a  re-examination  is  necessary. 

To  find  the  velocity  distribution,  an  additional  equation  is  required  since  there  are 
three  unknowns,  n,  v,  and  p,  appearing  in  Eq.  (12)  and  (13).  The  integral  continuity 
equation  is  used. 


-Li!  (23) 

Jr2  L 

Because  the  momentum  equation,  Eq.  (13),  is  nonlinear  in  u,  it  is  difficult  to  solve 
by  analytical  methods.  Instead,  a  direct  numerical  method  is  used  to  solve  this 
problem. 
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To  find  the  temperature  distributions  in  the  flowing  fluid  and  the  heat  fluxes 
between  the  flowing  fluid  and  the  PCM,  the  energy  equation  must  be  solved.  If  the 
velocities  and  the  boundary  conditions  for  the  energy  equation  are  known,  it  is  not 
difficult  to  solve  the  problem,  however,  if  the  boundary  condition  between  the  flowing 
fluid  and  the  PCM  is  unknown,  a  technique  of  finding  the  boundary  condition  is 
needed. 
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IV.  NUMERICAL  PROCEDURE 


A.  VELOCITY  PROFILES  IN  THE  ENTRANCE  REGION 


Although  some  methods  [17,18,19]  were  previously  used  to  find  the  laminar 
velocity  distribution  of  an  annular  passage  in  the  entrance  region,  the  previous 
methods  are  complicated  and  therefore  are  not  very  practical  for  this  study.  A  much 
simpler  method  will  be  used. 

By  choosing  an  arbitrary  value  for  pj,  then  have  only  two  unknowns,  u  and  v.  The 
procedure  of  finding  u  and  v  is:  (1)  give  an  arbitrary  value  for  v  (2)  solve  u  by  the 
Gaussian  elimination  method,  (3)  solve  v  by  the  Gaussian  elimination  method,  (4) 
compare  the  new  u  with  that  in  the  previous  iteration  and  if  the  difference  is  sufficiently 
small,  then  u  and  v  obtained  can  be  considered  satisfy  for  the  value  of  p.x,  otherwise, 
using  the  new  value  of  u  and  v  and  return  to  procedure  (2). 

Since  u  and  v  are  obtained  by  choosing  an  arbitrary  value  for  pj,  therefore  u  may 
or  may  not  satisfy  the  integral  continuity  equation.  An  iterative  solution  for  pj  is 
required  and  the  value  of  pj  is  considered  to  be  converged  if  and  only  if 


ru  dr 


l  -'2 


<  e 


(24) 


where  e  is  the  convergence  criteria.  The  strategy  to  find  pj  is  that  if  the  chosen  value 
of  pj  is  smaller  than  the  exact  value,  then 


r 


ru  dr  < 


1  -  r. 
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and  a  larger  value  of  pj  is  used  in  next  iteration.  The  Newton- Raphson  method  is  used 
to  find  pj. 

B.  ITERATIVE  SOLUTION  OF  THE  INNER-WALL  TUBE  TEMPERATURE 

Since  the  enthalpy  equation  nor  the  energy  equation  can  be  used  without  giving 
the  temperature  distribution  between  the  heat  transfer  fluid  and  the  PCM,  the  iterative 
solution  of  the  inner  tube  wall  temperature,  Tx,  becomes  a  starting  procedure  for  the 
calculation  of  the  heat  transfer  in  the  latent  thermal  storage  system.  Since  this 
temperature  is  between  the  inlet  fluid  temperature  and  the  initial  PCM  temperature, 
the  iterative  solution  converges. 


C.  SOLVING  THE  MELTING/FREEZING  PROBLEM 
The  boundary  conditions  of  Eq.(9)  are  given  by: 

,  PCp{Tx-Tf ) 


(25) 


d<}> 

dr 


=  0 

r= 0 


(26) 


d<f>  i 
dX  l*=o 


d<t> 
dX  'x=c 


=  0 


(27) 


Equation  (27)  is  based  on  the  assumption  that  the  two  end-plates  for  the  PCM  tube 
have  very  low  thermal  conductivities  so  that  they  are  considered  to  be  adiabatic. 
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The  enthalpy  distribution  in  the  PCM  at  a  fixed  axial  location  can  be  easily 
obtained  by  using  Eq.  (10)  and  the  appropriate  boundary  conditions.  The  temperature 
distribution  is  obtained  from  Eq.  (11).  The  new  0  from  the  last  time  step  is  a  result 
of  radial  conduction  and  axial  conduction  during  the  time  step.  The  0  obtained  from 
Eq.  (10)  is  not  appropriate  to  be  used  for  the  calculation  of  the  heat  flux  between  the 
heat  transfer  fluid  and  the  PCM.  The  equation  which  resembles  Eq.  (10)  but  lacks  the 
axial  conduction  term  is  given  below  and  is  used  for  the  purpose  of  calculating  the  heat 
flow. 


C1  -  05  2  K  Ku  -  *3 )  -  ( <t>?J  -  C.; )  r, 1 


At  ri+ ,  +  rt 


A  r 


(28) 


where 


0  = 


P(' 


_ 2 _ f 

Ax'  Jo 


p - - - dv 


The  transient  local  heat  flow  rate  for  a  length  of  Ajc'  at  r  =  rx  is  given  by 


YviPp\{9?+x  -  9?) 


A  t' 


(29) 


where 


Vi=n{r'2i+x-r'])  r?  Ax' 


and 
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A/'  = 


At  / 


2 

1 


Therefore  the  transient  local  heat  flow  per  unit  length  is  given  by 


Qxi(1’ x ') 


Ajt' 


(30) 


I).  TRANSIENT  TEMPERATURE  DISTRIBUTION  OF  THE  HEAT  TRANSFER 
FLUID 

From  the  integration  of  the  one-dimensional  heat  conduction  equation  for 
cylindrical  coordinates  with  the  boundary  conditions  where 


T  =  r,  at  r'  =  r\ 


and 


r  =  r2  at  r'  =  r'2 


and  using  Fourier's  law,  one  has 


Q  = 


2*M  r2  -  r,) 

In  (r'2lr\) 


(31) 


where  Q  is  the  heat  flow  rate  through  the  PCM  tube  per  unit  length. 

Equation  (31)  is  obtained  based  on  the  assumptions  that  the  tube  wall  is  very  thin 
so  that  the  axial  conduction  and  the  sensible  heat  absorbed  by  the  tube  are  neglected. 
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Since  the  heat  flow  through  the  PCM  tube  is  equal  to  the  heat  flow  from  the  PCM  tube 
to  the  PCM,  from  L'qs.  (30)  and  (31)  one  has 


r(!  &,(/.*')  ■litfrVr',) 

r2(/,x) - — - +  r  ,(,.*) 


(32) 


Upon  finding  the  temperature  distribution  of  the  HTF,  the  Nusselt  numbers  can 
then  be  calculated.  According  to  Fourier's  law,  the  local  heat  flow  per  unit  length 
between  the  heat  transfer  fluid  and  the  PCM  can  be  given  by 


Qxii1'  x')  =  2*  r'i  kh 


dT 

3r' 


r‘=  f  2 


(33) 


The  bulk  temperature  of  the  heat  transfer  fluid  is 


n- 


uTrdr 


ur  dr 


(34) 


By  definition,  the  Nussclt  number  is 


Nu  = 


hDh 

kH 


(35) 


The  heat  fiux  by  conduction  at  the  outer  PCM  tube  surface  is  equal  to  the  heat  flux 
by  convection  to  the  HTF.  Thus, 


h{Tb 


Ti  )  =  *» 


dr 


dr'  W=r’2 


The  combination  of  Eq.  (40)  and  Eq.  (41)  yields 


(36) 
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Nu  = 


(37) 


7(  »  ~  r2  )  37  j 

7),  —  Ti  dr  i  r=  r2 

It  is  remembered  that  all  the  values  concerning  the  heat  transfer  between  the  heat 
transfer  fluid  and  the  PCM  and  the  values  concerning  the  melting/freezing  in  the  PCM 
are  based  on  the  iterative  T,  and  it  may  not  be  converged.  Tx  is  considered  to  be 
converged  if  and  only  if 


I  Qx\  ~  Qx2  I  ^  6  (38) 

where  e  is  the  convergence  criteria.  If  T,  is  not  converged,  then  one  has  to  correct  it 
and  return  to  the  solution  procedure  described  of  solving  the  melting/freezing  problem. 
This  is  a  root  finding  problem  and  Tx  is  found  by  the  Newton- Raphson  method. 
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V.  EXPERIMENT 


The  experimental  apparatus  employed  in  this  study  is  shown  schematically  in  Fig. 
3.  The  experimental  setup  consists  of  three  systems:  the  instrumented  heat  exchanger 
system,  the  fluid  flow  control  system,  and  the  data  acquisition  system.  The 
instrumented  heat  exchanger  system  consists  of  a  video  camera,  a  paddlewhcel  flow 
meter,  controllable  light  source,  and  a  latent  thermal  storage  unit.  The  video  camera 
is  accompanied  by  a  video  cassette  recorder  (VCR)  for  photographic  data  recording. 
Specifications  for  the  video  camera  arc  provided  in  Table  1.  The  paddlewhcel  flosensor 
is  accompanied  by  a  digital  meter  for  heat  transfer  fluid  flow  rate  monitoring. 

Table  1.  VIDEO  CAMERA  SPECIFICATIONS 

Unit  Panasonic  WV-3250/8AF  color  video 

camera 

Specification  (lens) 

Auto  focus 

Zoom  ratio  8x 

Focal  length  10.5-84  mm 

M ax.  apeture  ratio  1 : 1 .4  ( F 1 .4) 

Macro  range  5  cm  -  1.3  m 

White  balance  Auto  tracing  white  balance  or 

manual  setting 

The  shell  and  tube  heat  exchanger  consists  of  a  3  in.  (7.62  cm)  inner  diameter 
shell  around  a  removeable  cartridge  of  tubes.  Each  cartridge  consists  of  eight  capped 
tubes  filled  with  PCM  and  two  support  disks.  The  disks  were  designed  to  provide  four 
functions,  which  are:  provide  a  support  frame  for  the  tube  bundle,  align  fluid  flow  over 
the  tubes,  align  thermocouples  and  form  a  plenum  region  between  the  cartridge  and  the 
end  plates  of  the  shell.  The  plenum  regions  serve  to  maintain  a  full  level  of  fluid  flow 
over  the  PCM  containment  tubes.  Both  the  shell  and  cartridge  were  fabricated  from 
Lexan.  Lexan  is  a  transparent  polycarbonate,  thereby  enabling  visual  access  to  the 
phase  change  phenomena.  As  a  heat  transfer  material  Lexan  suffers  from  low  thermal 
conducitvity,  but  the  limiting  heat  transfer  resistance  in  this  shell  and  tube  heat 
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Figure  3.  Experimental  Apparatus 


exchanger  was  determined  to  be  the  PCM.  Therefore  for  experimental  purposes  Lcxan 
was  an  excellent  choice.  The  heat  exchanger  is  sealed  with  end  plates  that  are  grooved 
for  a  gasket  to  prevent  leakage.  These  end  plates  have  holes  tapped  for  the  inlet  and 
outlet  fluid  flow  ports,  and  also  for  the  thermocouple  fittings.  A  total  of  four 
copper-constantan  subminiature  thermocouples  were  fitted  in  the  inlet  and  outlet 
positions  of  the  cartridge.  Another  copper-constantan  thermocouple  was  fabricated 
and  placed  inside  one  of  the  PCM  containment  tubes  approximately  halfway  between 
the  inlet  and  outlet  for  PCM  temperature  monitoring. 

A  second  shell  and  tube  heat  exchanger,  utilizing  the  same  shell  from  the  first 
heat  exchanger,  was  designed  to  improve  photographic  capabilities.  A  single  2  in.  (5.08 
cm)  diameter  PCM  containment  tube  was  fabricated  out  of  Lexan,  with  a 
thermocouple  fitting  located  at  the  front  end  of  the  tube.  A  copper-constantan 
thermocouple  was  used  to  monitor  the  PCM  temperature,  at  an  axial  distance  from  the 
inlet  of  the  PCM  tube  of  8.4  in.  (21.34  cm).  Incorporated  into  the  design  was  a  higher 
outlet  fluid  flow  port  to  facilitate  the  purging  of  the  initial  air  inside  the  storage  unit. 
Table  2  lists  the  specifications  of  both  heat  exchanger  configurations;  the  tube  bundle 
and  the  single  containment  tube. 

Table  2.  SHELL  AND  TUBE  HEAT  EXCHANGER  SPECIFICATIONS 


Parameter 

Tube  Bundle 

Single  Tube 

Shell  Inner  Diameter  (cm) 

7.62 

7.62 

Shell  Wall  Thickness  (cm) 

0.47 

0.47 

Tube  Outer  Diameter  (cm) 

0.953 

5.08 

Tube  Inner  Diameter  (cm) 

0.635 

4.76 

Material 

Lexan 

Lexan 

Conductivity  (W/mC) 

7.96 

7.96 

Total  PCM  Volume  (cm3) 

1  5.6 

882.0 

Flow  Length  (cm) 

49.2 

49.5 
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The  fluid  (low  control  system  consists  of  a  fin  and  tube  air  cooled  heat  exchanger 
(cooler),  a  constant  temperature  both  circulator,  an  electronic  temperature 
programmer,  and  a  remote  sensor.  The  bath/circulator  incorporates  a  dual  pump,  a 
heater,  and  a  sensor  into  a  fluid  reservoir.  An  electronic  controller  is  mounted  on  top 
of  the  reservoir.  The  cooler  is  used  in  conjunction  with  the  bath/circulator  during  the 
freezing  process.  The  cooler  is  also  used  in  the  temperature  programmer  is  a  linear 
programmmer  which  is  connected  to  the  bath/circulator  via  the  remote  sensor.  The 
programmer  can  t  "uratcly  control  the  temperature  of  the  heat  transfer  fluid.  Table 
3  gives  the  details  of  the  fluid  flow  control  system. 


Table  3.  FLUID  FLOW  CONTROL  SYSTEM  SPECIFICATIONS 

Constant  Temperature  Bath  Circulator: 

Unit  Neslab  EX-100DD 

Specifications 

Range  Ambient  to  266  F  ( 1 30  C) 

Stability  0.018  F  (0.01  C) 

Max  Pumping  3.35  gpm  (13  1/min) 

Heater  800  Watts 


Fin  and  Tube  Air  Cooler  Heat  Exchanger  (Cooler): 
Unit  Neslab  Aerocool  100 

Specifications 

Lower  limit  3.6  F  (2  C)  above  ambient 

Fan  speed  adjustable 


Remote  Sensor: 
Unit 

Specifications 

Sensor 

Resolution 


Neslab  RS-1 

Platinum  RTD 
0.18  F  (0.1  C) 


Electronic  Temperature  Programmer: 

Unit  Neslab  ETP-3 

Specifications 

Programming  Linear 

Time  Interval  10-10,000  minutes 

Clock  Pulse  100  microseconds 


The  data  acquisition  system  is  a  microcomputer  based  data  logger.  This  system 
consists  of  an  IBM  PC/XT  with  a  Strawberry  Tree  Analog  Connection  System 
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installed,  a  high  resolution  color  monitor,  a  high  quality  (HQ)  VCR,  and  a  digital  fluid 
flowmeter. 

The  Strawberry  Tree  Analog  Connection  has  a  menu  driven  data  acquisition 
program  which  also  has  process  control  capabilities.  Connections  for  the 
thermocouples  and  for  the  signal  conditioner  of  the  flowmeter  are  provided  by  an 
analog  to  digital  (A/D)  terminal  box.  The  terminal  box  provides  input  channels  to  the 
microcomputer  for  on-screen  flow  rate  and  temperature  monitoring  and  floppy  disk 
recording.  The  monitor  has  dual  capability  as  a  high  resolution,  digital  RGB  input 
monitor  and  as  a  color  display  for  the  VCR.  The  HQ  VCR  is  used  to  monitor,  record 
and  replay  the  photographic  image  of  the  solid-liquid  phase  front  motion.  The  VCR 
has  a  slow  speed  replay  feature  for  frame-by-frame  advance  of  the  type.  The  VCR 
recording  can  be  still  framed  and  the  resulting  picture  photographed  by  a  35mm 
camera.  The  digital  fluid  flowmeter  converts  the  frequency  signal  from  the  paddlewheel 
flosensor  to  a  four  digit  display  of  the  volumetric  flow  rate.  Table  4  gives  the  details 
of  the  data  acquisition  system. 


Table  4.  DATA  ACQUISITION  SYSTEM  SPECIFICATIONS 


Microcomputer  based  Data  Logger: 

Unit  IBM  PC/XT  with  Strawberry  Tree  Analog 

Connection  System 


Specifications 
No.  of  Channels 
Resolution 
Accuracy 
Read  rate 


8 

0.015% 

0.04% 

0.3  to  5  msec 


Video  Cassette  Recorder: 

Unit  Toshiba  M-2400  VHS  video  cassette 

recorder 

35mm  Camera  &  Lens: 

Camera  Minolta  Advanced  Maxxum  7000 

35  mm  SLR  camera 


Lens 


Minolta  Maxxum  70-210  mm  f/4.0 
with  macro  zoom 


24 


The  video  camera  and  VCR  were  employed  in  the  experiments  because  of 
unsuccessful  attempt  to  use  a  Micorneye  Bullet  digital  imaging  system  to  track  the 
phase  front  progression.  Attempts  were  made  to  improve  the  digital  image  that  was 
received  from  the  system.  Modifications  to  the  background  lighting,  focal  lenght  of  the 
lens,  and  the  apeture  of  the  lens  were  made.  However  these  modifications  could  not 
overcome  the  depth  of  field  limitaions  of  the  Microneye  Bullet  system  when  viewing  the 
curved  surfaces  of  the  PCM  containment  tubes.  The  depth  of  field  of  the  Microneye 
camera  was  less  than  1  cm,  which  was  too  small  to  obtain  an  accurate  image.  Table 
5  is  a  comparision  between  the  Microneye  digital  imaging  system  and  the  VCR/video 
camera  setup  used  in  the  experiments.  The  video  camera  provides  a  greater  depth  of 
field  and  faster  iris  control  when  compared  to  the  Microneye  camera,  and  this  led  to  a 
much  improved  image  of  the  PCM  phase  front. 


A.  TEST  MATERIALS 

The  phase  change  materials  were  selected  based  on  their  availability  and  melting 
temperature.  Low  melting  temperature  materials  were  selected  because  of  their 
application  in  Spacecragr  Thermal  Management  Systems.  Importance  was  also  placed 
on  materials  that  would  allow  clear  visual  indentification  of  the  phase  change  front 
motion. 

PI  16,  a  Sun  wax,  was  selected  as  a  PCM  to  be  tested.  PI  16  is  a  paraffin  which 
is  clear  as  a  liquid  but  opaque  as  a  solid.  Melting  occurs  with  limited  segregation  of 
components  and  upon  freezing  PI  16  tends  to  pull  away  from  it's  container,  which  has 
an  apperciable  effect  on  heat  transfer.  Also,  paraffin  waxes  display  little  supercooling. 

E-IOOO,  a  polyethylene  glycol,  is  a  stable  organic  polymer.  E-1000  displays  no 
supercooling  and  is  compatible  with  many  materials.  E-1000  is  white  as  a  solid  and 
transparent  as  a  liquid. 
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C-31  is  a  cross-linked  polymer  formulation  of  Na2S04  —  10H2O.  The  phase 
change  process  is  not  the  traditional  freeze/melt  concept  but  a  solid  transition. 
Therefore,  there  is  no  phase  change  driven  separation  or  stratification.  Properties  for 
the  test  materials  may  be  found  in  Table  6. 

Water  was  selected  as  the  heat  transfer  fluid.  It  is  clear,  abundantly  available, 
and  inexpensive.  It  is  used  in  a  variety  of  thermal  storage  and  heat  transfer  related 
applications.  Also  water  is  transparent  and  thereby  does  not  interfere  with  the  viewing 
of  the  phase  change  phenomena. 
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Table  5.  COMPARISON  OF  IMAGING  TECHNIQUES 
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Table  6.  PROPERTIES  OF  THE  TEST  MATERIALS 


PROPERTIES 

PCM 

P'1  1  6 

C  — 31 

E— 1 000 

Melting  Point  (C) 

43.0 

31  .0 

38.0 

Heat  of  Fusion  (KJ/Kg) 

226 

200 

1  72 

Specific  Heat  (KJ/KgC) 

2.75 

2.09 

2.05 

Conductivity  (W/mC) 

0.24 

0.61 

Density— Solid  (Kg/m3) 

818 

1  420 

1  123 

Density— Liquid  (Kg/m3) 

760 
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VI.  EXPERIMENTS 


A.  DIFFERENTIAL  SCANNING  CALORIMETRY 


Differential  scanning  Calorimetry  (DSC)  tests  were  performed  on  pi  16  to  verify 
melting  and  recrystallization  temperatures,  and  the  latent  heat  of  fusion  value.  A 
Pcrkin-EImer  DSC-4  Differential  Scanning  Calorimeter  and  a  Perkin-EImer  Model 
3600  Data  Station  erer  employed  to  perform  the  tests.  Figures  4  and  5  represent  the 
results  of  these  tests. 

DSC  tests  were  also  performed  on  C-31.  However,  difficulties  were  encountered. 
C-3I  is  a  polymer  formulation  of  a  salt  hydrate  and  during  the  tests  the  water  would 
be  driven  off,  therefore  giving  ficitious  results  of  the  test  sample  energy  output. 
Mechanical  sealing  ot  the  sample  containers  was  ineffective  in  preventing  this 
phenomena.  Chemical  bonding  of  the  sample  containers  was  also  ineffective. 

B.  PI  16  SUN  WAX 


PI  16,  a  Sun  wax,  was  tested  using  two  different  heat  exchanger  assemblies.  A 
cartridge  design  and  a  single  containment  tube  exchanger  was  used. 

1.  Tube  Bundle  Results.  Experiments  were  performed  for  both  the  melting  and 
freezing  of  the  PCM.  In  each  case  two  trials  were  completed  at  each  of  two  flow  rates: 
0.6  gpm  (2.27  1/min)  and  1.22  gpm  (4.62  1/min). 

For  the  melting  process,  data  acquisition  began  once  thermal  equilibrium  in  the 
heat  exchanger  was  cxtablished  at  a  temperature  below  the  melting  point  of  the  PCM. 
The  setpoint  of  the  constant  temperature  bath  wai  set  at  48  C  and  inlet,  outlet  and 
PCM  temperature  were  recorded  until  melting  was  complete.  For  the  freezing  process, 
data  acquisition  began  once  the  setpoint  of  the  constant  temperature  bath  was  set  at 
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MCAL/SEC  „  MCAL/SEC 


% 


Figure  5. 


raaw  rae.wAxa.D4  temperature  <c>  DSC 

DATE.  87/07/22  TIME.  10. 35 


none  rae.xAX3.04  temperature  cc>  DSC 

OATE.  87/07/20  TIME)  22.49 


DSC  Results  for  Melting  (A)  and  Freezing  (B)  of  PI  16 
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43  C,  the  air  cooler  was  turned  on,  and  data  acquisition  was  continued  until  freezing 
was  complete. 

Figures  6-9  represent  on  selected  trial  for  each  experimental  condition.  The  inlet 
temperature  repersents  the  temperature  of  a  single  Thermocouple  located  at  the 
beginning  of  the  PCM  tubes.  The  PCM  temperature  represents  the  temperature  of  one 
of  the  two  thercocouples  located  interstitially  in  the  PCM.  The  outlet  temperature  is 
the  temperature  of  one  of  the  three  thermocouples  located  at  the  end  of  the  PCM 
tubes.  This  method  of  data  reduction  was  chosen  because  of  the  lack  of  data  scatter 
between  trials  and  also  due  to  the  lack  of  data  scatter  between  thermocouples  in  the 
outlet  and  PCM  for  each  trial.  As  a  measure  of  repeatability,  heating  and  cooling  rates 
of  the  PCM  for  each  trial  were  calculated.  F  or  the  low  flow  rate,  heating  rates  of  0.28 
C/min  and  0.74  C/min,  and  cooling  rates  of  0.44  C/min  were  found.  Comparable 
results  were  obtained  for  the  high  flow  rate.  Reynolds  and  Stefan  numbers  were 
calculated  in  all  experiments  as  follows: 

Re  =  V  Dh  p  I  u  and  Ste  =  cp  AT /  h 

where  v  is  the  velocity  of  the  heat  transfer  fluid,  Dh  is  the  hydraulic  diameter,  p  is  the 
density  of  the  HTF,  p  is  the  viscosity  of  the  HTF,  cp  is  the  average  specific  heat  of  the 
PCM,  AT  is  the  temperature  swing  of  the  PCM,  and  h  is  the  latent  heat  of  fusion.  The 
temperature  swing  of  the  PCM  in  this  case  would  be  the  difference  of  the  final  PCM 
temperature  and  the  initial  PCM  temperature. 
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FREEZING  OF  PII61  R 


Figure  8.  Experimental  Temperature  for  Freezing  of  PI  16,  0.6  gpm  (2.27  1/min) 
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Figure  9.  Experimental  Temperature  for  Freezing  of  PI  16,  1.22  gpm  (4.62  1/min) 


Figures  6  and  7  represent  selected  results  for  the  melting  process  for  the  low  and 
high  flow  rates  respectively.  The  PCM  temperature  was  not  constant  during  melting. 
This  is  most  likely  due  to  the  non-uniformity  of  PI  16.  The  inlet  temperature  is  shown 
to  overshoot  the  constant  temperature  bath  setpoint  of  48  C,  as  the  bath  compensates 
for  the  thermal  inertial  of  the  heat  exchanger.  Figures  8  and  9  represent  selected  results 
for  the  freezing  process  for  the  low  and  high  flow  rates  respectively.  The  inlet 
temperature  decreases  linearly  without  overshooting  the  bath  setpoint  of  43  C.  The 
heater  in  the  constant  temperature  bath  and  the  air  cooler  work  together  to  decrease 
the  temperature  without  overshooting  the  setpoint.  It  is  important  to  note  that  the 
air  cooler  is  not  a  refrigeration  unit  and  that  the  maximum  amount  of  cooling  was 
provided  during  the  freezing  experiments. 

The  video  camera  provided  a  qualitative  insight  to  the  nature  of  melting  and 
freezing  of  PI  16.  PI  16  wax  is  clear  as  a  liquid  but  opaque  as  a  solid,  thus  enabling 
visual  identification  of  the  phase  front.  For  each  trial,  the  video  camera  and  VCR  were 
used  to  tape  the  entire  experiment.  At  the  low  flow  rate,  visulaization  was  limited  dure 
to  air  bubbles  adfering  to  the  outside  of  the  PCM  tubes.  Air  bubbles  were  not  a 
problem  at  the  high  flow  rate  as  melting  progressed  non-uniformly  in  the  video. 
Observation  of  the  freezing  process,  however,  was  limited  due  to  the  initial  frosting  of 
the  inside  PCM  tube  walls  at  the  onset  of  freezing. 

2.  Single  Tube  Results.  An  experiment  was  performed  for  the  melting  of  PI  16 
in  a  single  containment  tube  at  a  flow  rate  of  0.65  gpm  (2.95  1/m).  Thermal  equilibrium 
was  established  in  the  heat  exchanger  at  a  temperature  of  39  C.  Flow  was  shut  off  to 
the  heat  exchanger.  Once  the  bath  temperature  reached  49  C,  fluid  flow  was  initiated 
into  the  heat  storage  unit  and  the  data  acquisition  process  was  begun  and  continued 
unitl  melting  was  complete.  During  the  experiment,  the  inlet  temperature  was  held 
constant  at  49  C.  Figure  10  represents  one  test  for  the  melting  of  PI  16  in  the  single 


tube.  The  inlet,  outlet  and  PCM  temperatures  are  that  of  single  thermocouples.  A 
discontinuity  in  time  of  the  PCM  temperature  is  evident.  This  represents  the  phase 
change  transition.  Figure  11  and  12  are  a  succession  of  photographs  of  the  phase 
change  process  for  PI  16  was  in  a  single  tube.  In  each  photograph  the  time  is  displayed. 
The  thermocouple  is  visible  and  it  acts  to  hold  the  solid  PCM  from  falling  to  the 
bottom  of  the  tube.  The  other  photographs  display  the  iowerhalf  of  the  PCM  tube. 
The  time  displayed  on  the  photographs  is  consistent  with  the  time  displayed  on  the 
x-axis  of  Fig.  10.  Additional  tests  for  P 1 16  were  not  performed  due  to  the  tube  failure 


Figure  10.  Experimental  Temperature  for  Melting  of  PI  16,  0.65  gpm  (2.46  1/min) 
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Figure  1 1.  Experimental  Photographs  for  Melting  of  PI  16,  0.65  gpm  (2.46  lAnin) 


Figure  12. 


C.  CALOTHERM  C-31 


Experiments  were  performed  for  the  melting  of  C-31  in  the  single  containment 
tube  at  a  flow  rate  of  0.3  gpm  (1.34  1/min).  Thermal  equilibrium  was  extablished  in  the 
heat  exchanger  at  a  temperature  of  30  C.  Fluid  flow  was  shut  off  to  the  heat 
exchanger.  Once  the  bath  temperature  reached  55  C,  fluid  flow  was  initiated  into  the 
heat  exchanger  and  the  data  acquisition  process  was  begun  and  continued  until  melting 
was  complete.  When  the  experiment  was  in  progress,  the  inlet  temperature  was  held 
constant  at  55  C. 

Figure  13  displays  two  selected  trials  of  the  experiments  performed  on  C-31.  The 
inlet,  outlet  and  PCM  temperatures  represent  the  temperature  of  single  thermocouples. 
The  PCM  temperature  increase  gradually,  and  no  sharp  melting  point  transition  is 
noticed.  This  is  a  result  of  the  fact  that  phase  change  is  actually  a  solid-solid  transition. 
The  top  photograph  in  Fig.  14  displays  stress  cracks  in  the  Lexan  of  the  single 
containment  tube  that  were  presummably  induced  by  the  volume  change  of  C-31.  The 
damaging  volume  change  was  not  expected  since  the  phase  change  is  a  solid  transition. 

D.  E-1000 


Attempts  were  made  to  experiment  with  E-1000  contained  in  the  eight  cylindrical 
tubes  of  the  cartridge  design.  Upon  one  melting-freezing  cycle,  tube  failure  occurred. 
Besides  the  leaking  of  was,  volume  change  with  phase  change  induced  stress  was 
evident.  The  bottom  photograph  in  Fig.  18  displays  the  stress  cracks. 
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Figure  14.  Photographs  of  Stress  Cracks  on  Lexan  Tubes 
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VII.  VERIFICATION  OF  ACCURACY 


Since  in  recent  studies  many  of  the  methods  used  are  cither  new  or  improved 
from  other  methods,  a  comparison  of  the  recent  numerical  results  with  experimental 
results  or  the  results  obtained  by  others  is  very  important.  This  not  only  avoids 
mistakes  made  in  recent  studies  but  also  provides  an  examination  of  the  accuracy  of 
recent  methods. 

We  found  that  many  of  the  previous  studies  were  restricted  to  limited  boundary 
conditions  or  simpler  cases,  therefore  the  comparisons  with  previous  works  are  quite 
limited. 

A.  MELTING/FREEZING  OF  PCM  IN  A  HORIZONTAL  TUBE. 

A  comparison  of  the  experimental  data  of  Viskanta  [  22  ]  and  the  predictions 
by  the  enthalpy  method  with  n-heptadecane  as  the  PCM,  given  in  Fig.  15,  shows  quite 
good  agreement.  In  Fig.  15  the  Stefan  number,  Ste,  and  the  initial  dimensionless 
temperature,  <f>ini4  are  defined  as: 


Ste 


cP{Tf- r,) 
A 


(39) 


r-r, 

9m  77-  r, 


(40) 


B.  ENTRANCE  REGION  VELOCITY  DISTRIBUTION 

The  method  used  in  this  investigation  to  find  the  laminar  velocity  distribution  in 
the  entrance  region  of  an  annulus  is  very  simple  when  compared  with  the  eigenvalue 
method.  When  computing  on  IBM  4381  system  for  41  x  501  grids,  the  CPU  time  is 
only  86  seconds.  The  comparison  of  maximum  velocity  distribution  by  1  Icato'  ct  al. 
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Figure  15.  Comparison  of  predicted  solid-liquid  interface  position  with  data  for  solidification  of  n-heptadecanc 


[19]  and  this  work  is  given  in  Fig.  16.  When  the  flow  is  fully  developed,  the  velocity 
profile  obtained  by  this  method  has  a  maximum  difference  of  only  0.02%  compared 
with  Eq.  (22)  for  the  case  of  r2  =  0.5  and  1 1  grid  points  in  the  radial  direction.  A 
comparison  of  the  fully  developed  velocity  profiles  obtained  from  Eq.  (22)  and  from  the 
numerical  method  for  r2  =  0.5  is  shown  in  Table  7. 

C.  HEAT  TRANSFER  IN  ANNULAR  PASSAGES 

Comparisons  of  the  results  for  heat  transfer  in  annular  passages  for  flow  in  the 
fully  developed  region  with  Lundberg's  [13]  and  Worsoe-Schmidt's  [  14  ]  are  given  in 
Table  8.  For  flow  in  the  entrance  region,  the  results  compared  with  Heaton's  [  19  ]  for 
the  case  of  uniform  heat  flux  at  the  inner  wall  are  shown  in  Table  9.  In  Heaton's 
analysis,  the  effects  of  v  on  the  flow  field  are  underestimated  because  of  linearization, 
therefore  the  Nu  numbers  obtained  by  Heaton  do  not  agree  very  well  with  the  present 
results.  A  better  agreement  is  obtained  by  neglecting  v  of  Eq.  (14). 
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Figure  16.  Maximum  velocity  variation  in  the  entrance  region  of  an  annulus 


Table  7.  VELOCITY  DISTRIBUTIONS  IN  FULLY  DEVELOPED  REGION 
FOR  r2  =  0.5 


r-r2 

Numerical  method 

11  grids  41  grids 

Equation  (22) 

1  —  r2 

0.1 

0.60266 

0.60278 

0.60279 

0.2 

1.03899 

1.03912 

1.03913 

0.3 

1.32605 

1.32616 

1.32616 

0.4 

1.47716 

1.47721 

1.47721 

0.5 

1.50284 

1.50283 

1.50283 

0.6 

1.41159 

1.41154 

1.41153 

0.7 

1.21036 

1.21029 

1.21028 

0.8 

0.90493 

0.90485 

0.90485 

0.9 

0.50011 

0.50006 

0.50006 
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Tabic  8.  SOLUTIONS  OF’  I  LOW  IN  FULLY  DUVHLOPHD  RFGION  FOR 


r2  =  0.5 


X  NU 

Present  Lundberg  Worsoe-Schmidt 


Constant  Heat  Flux  at  Inner  Wall 


0.0001 

34.985 

34.6 

34.233 

0.0010 

16.39 

16.4 

16.356 

0.0025 

12.378 

12.37 

12.367 

0.005 

10.132 

10.13 

10.127 

0.010 

8.435 

8.433 

.... 

0.025 

6.932 

6.931 

— - 

0.05 

6.353 

6.353 

— - 

0.1 

6.191 

6.192 

.... 

Constant  Temperature 

at  Inner  Wall 

0.0001 

29.802 

28.4 

28.456 

0.0010 

13.760 

13.68 

13.701 

0.3025 

10.445 

10.42 

10.427 

0.005 

8.610 

8.602 

8.603 

0.010 

7.249 

7.246 

7.247 

0.025 

6.118 

6.117 

0.05 

5.784 

5.785 

— - 

0.1 

5.738 

5.738 
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Table  9.  SOLUTIONS  OF  UNIFORM  HEAT  FLUX  AT  INNER  WALL  FOR 


FLOW  IN  ENTRY  REGION  WITH  Pr  =  10 


X 

Nu 

Present 

Neglect  v 

Heaton 

0.0010 

16.081 

16.468 

16.68 

0.0025 

12.255 

12.535 

12.60 

0.005 

10.079 

10.187 

10.20 

0.01 

8.414 

8.453 

8.43 
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VIII.  RESULTS  AND  DISCUSSION 


Figure  17  shows  the  Nusselt  numbers  in  the  thermal  entrance  region  for  two 
cases,  one  by  considering  that  the  flow  is  fully  developed  and  another  case  by  assuming 
that  the  flow  field  is  developing.  We  see  that  near  the  inlet,  Nusselt  numbers  of  the 
developing  flow  field  are  larger  than  those  of  the  fully  developed  flow  field.  This  is  due 
to  the  fact  that  in  the  velocity  developing  region,  the  boundary  layer  is  thinner, 
therefore  the  drag  force  is  larger  and  this  makes  the  Nusselt  numbers  larger.  As  the 
HTF  flows  down  from  the  inlet,  the  Nusselt  numbers  of  the  two  different  flow  fields 
become  closer. 

From  Eq.  (13)  and  Eq.  (14),  one  can  conclude  that  the  Reynolds  number  will  not 
affect  the  local  heat  flow  rates.  This  is  effective  only  when  X  is  defined  by 
X—  x '  I  Re  Pr  Dh.  The  effects  of  changeing  the  initial  temperature  of  the  HTF  and  the 
PCM  and  changing  the  inlet  temperature  of  the  HTF  of  the  thermal  storage  unit  on 
the  Nusselt  number  are  shown  to  be  insignificant. 

Figure  18  shows  the  transient  temperature  distribution  of  the  PCM  for  flow  in 
the  entrance  region.  The  initial  temperature  of  the  PCM  and  HTF  is  38°C.  The  inlet 
temperature  of  the  HTF  is  53° C.  Since  the  temperature  of  the  PCM  is  initially  lower 
than  the  melting  point  temperature,  the  solid  PCM  temperature  increases  to  the 
melting  point  and  melting  occurs.  However,  as  time  increases,  the  solid  PCM  will 
approach  the  melting  temperature,  a  constant  temperature.  The  temperature  gradient 
of  the  PCM  at  the  inner  tube  wall  will  decrease  with  time.  This  means  that  the  heat 
flow  rate  between  the  HTF  and  PCM  will  decrease  with  time. 

The  transient  temperature  distributions  of  the  HTF  are  shown  in  Figs.  19  and 
20.  Near  the  inlet,  the  temperature  of  the  HTF  is  uniform  except  at  the  region  which 
is  near  the  outer  wall  of  the  PCM  tube.  For  the  downstream  region,  the  temperature 
of  the  HTF  at  the  center  portion  between  r2  and  rQ  is  higher  than  that  near  the  tube 
wall.  This  phenomenon  is  due  to  the  parabolic  velocity  distribution.  However,  as  time 
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Figure  17.  EfTect  of  developing  velocity  profiles  on  Nussclt  number 


.  15 


<VJ 

at 

CD 

cn 

o 

m 

CO 

rl 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

l  l 

<t> 


51 


Figure  18.  Temperature  d 
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increases,  the  temperature  gradient  appears  to  be  nonzero  only  at  the  region  which  is 
near  to  r2.  From  Figs.  19  and  20  one  would  discover  that  as  X  decreases,  the 
temperature  gradient  at  r2  increases.  Thus,  largei  heat  flow  rates  are  predicted  in  the 
entrance  region  than  in  the  fully  developed  region  downstream.  Note  that  when  time 
increases,  the  local  heat  flow  rates  decrease. 

For  melting  of  the  PI  16  wax  in  the  latent  heat  storage  unit,  both  experimental 
data  and  predicted  results  showed  fair  agreement.  However,  after  the  experiment 
progressed  100  minutes,  air  bubbles  generated,  which  had  significant  influence  on  the 
experimental  results,  and  therefore  comparison  between  the  experimental  data  and 
predicted  results  becomes  difficult.  During  the  initial  melting  period,  a  minimal 
amount  of  slumping  of  the  solid  PCM  was  boserved.  The  end  plugs  of  the  PCM 
capsules  offered  sufficient  thermal  resistance  ot  allow'  the  solid  PCM  to  remain 
attached  to  these  plugs  for  a  substantial  period  of  the  melting  process. 

The  initial  time  evolution  of  the  solid-liquid  interface  positions  was  plotted  on 
Fig.  21  for  constant  Renolds  and  Prandtl  numbers  of  683  and  3.93,  respectively.  The 
figure  clearly  shows  that  in  the  upstream  region  of  the  HTF,  the  melting  of  the  PCM 
is  faster  than  in  the  downstream  region.  This  is  due  to  the  fact  that  in  the  upstream 
region,  the  temperature  of  the  HTF  is  higher  and  both  the  temperature  and  velocity 
profiles  are  developing.  Figure  21  also  indicates  that  the  melting  is  slower  as  time 
increases,  and  this  ch?rHcteri7es  that  the  heat  flow  rates  between  the  HTF  and  PCM 
decrease. 
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Figure  21.  Time  evolution  of  the  solid-liquid  interface  positions  for  melting  of  PI  16  wax 


IX.  CONCLUSIONS 


Analysis  and  experiments  have  been  performed  for  a  heat  transfer  study  of  a 
latent  heat  storage  unit.  Melting/freezing  problems  were  solved  by  a  finite-difference 
method  based  on  the  enthalpy  formulation.  Convection  problems  of  the  HTF  were 
solved  by  first  obtaining  the  velocity  distribution  and  then  by  solving  the  energy 
equation.  In  addition  to  the  continuity  equation  and  the  momentum  equation,  the 
integral  continuity  equation  was  also  used  for  calculating  the  velocity  distributions. 

From  the  obtained  results  one  can  conclude  that:  (1)  The  developing 
temperature  profiles  increase  the  Nusselt  number  significantly.  The  developing  velocity 
profiles  can  also  increase  the  Nusselt  number.  However,  the  increase  is  not  so 
significant  when  compared  with  the  increase  due  to  the  developing  temperature 
profiles.  (2)  Near  the  inlet  region,  the  effect  of  r2  on  the  Nusselt  number  can  not  be 
neglected.  The  Nusselt  number  increases  with  the  increase  of  r2.  In  the  far 
downstream  region,  the  effect  of  r2  on  the  Nusselt  number  is  small  and  can  be 
neglected. 
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LIST  OF  ACRONYMS 


A/D 

Analog  to  Digital 

DSC 

Differential  Scanning  Calorimetry 

HQ 

High  Quality 

HTF 

Heat  Transfer  Fluid 

PCM 

Phase  Change  Material 

SCTMS 

Spacecraft  Thermal  Management  Systems 

VCR 

Video  Cassette  Recorder 
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A, 
A  i, 

A2, 

a3, 

B, 
c, 

C, 
D', 
D'h, 
K 

i, 

k, 

L\ 

M, 

P'< 

Q, 

Qx  i- 
Qx  2- 

r', 

r, 

t, 

u, 

u\ 

V', 


area 

constant 

constant 

constant 

constant 

specific  heat 

constant 

diameter 

hydraulic  diameter,  2 (r'a  —  r'2) 

heat  transfer  coefficient 

specific  enthalpy 

thermal  conductivity 

length  of  tube 

constant 

pressure 

heat  flow  per  length 

local  heat  flow  per  length 

local  heat  flow  per  length 

heat  flow  rate 

spatial  coordinate 

temperature 

time 

specific  internal  energy 
axial  component  of  velocity 
radial  component  of  velocity 
spatial  coordinate 
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Greek  symbols 


a, 

thermal  difTusivity 
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latent  heat  of  fusion 

P, 

density 

o, 

kinematic  viscosity 

e  , 

convergence  criteria 

A, 

step  size 

Dimensionless  quantities 

Fo, 

Fourier  number,  dimensionless  time, 

L, 

length  of  tube,  L'  /  Re  Pr  D’0 

Nu, 

Nusselt  number 

P . 

pressure,  p'  /  (p  u'2m  /  2) 

Pr, 

Prandtl  number,  u/oc 

Re, 

Reynolds  number,  u‘m  D’h  /  o 

Re’, 

Reynolds  nu  nber,  2u'm  r'0  /  v 

r. 

spatial  coordinate,  r'jr'0 

Ste, 

Stefan  number,  Eq.(39) 

T, 

temperature,  (T  -  T j)!{T e-  T j) 

M, 

axial  component  of  velocity 

radial  component  of  velocity 

X, 

spatial  coordinate,  x’/Re  Pr  D'h 

x, 

spatial  coordinate,  x'/Re'  r'0 
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P. 

t, 

1, 

2, 


spatial  coordinate, 
spatial  coordinate,  x’/Re  r’0 
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scaled  dimensionless  temperature,  Eq.(8) 
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mean  bulk  value 
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fusion 

heat  transfer  fluid 
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phase  change  material 
tube 
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